ConvertGeodeticToHotine Subroutine

SUBROUTINE ConvertGeodeticToHotine &
!
(lon, lat, k, lon0, lat0, azimuth, a, e, falseN, falseE, x, y)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString


IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (IN) :: lat !!geodetic latitude [radians]
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: lon0 !!longitude of center [radians]
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of center [radians]
REAL (KIND = float), INTENT (IN) :: azimuth !!azimuth of centerline
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: falseN !!false northing
REAL (KIND = float), INTENT (IN) :: falseE !!false easting


!Arguments with intent (out):
REAL (KIND = float), INTENT (OUT) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (OUT) :: y !!northing coordinate [m]

!Local declarations:
REAL (KIND = float) :: e2, e4, e6, e8
REAL (KIND = float) :: B
REAL (KIND = float) :: B1
REAL (KIND = float) :: t0
REAL (KIND = float) :: D
REAL (KIND = float) :: D2
REAL (KIND = float) :: F
REAL (KIND = float) :: H
REAL (KIND = float) :: G
REAL (KIND = float) :: g0
REAL (KIND = float) :: l0
REAL (KIND = float) :: uc
REAL (KIND = float) :: gc
REAL (KIND = float) :: u
REAL (KIND = float) :: v
REAL (KIND = float) :: Q
REAL (KIND = float) :: S
REAL (KIND = float) :: TT
REAL (KIND = float) :: UU
REAL (KIND = float) :: VV
REAL (KIND = float) :: t
REAL (KIND = float) :: c

!------------end of declaration------------------------------------------------

!Eccentricities squared
e2 = e**2.
e4 = e2**2.
e6 = e**6.
e8 = e**8.

!calculate constants
B = ( 1. + e2 * (COS (lat0))**4. / (1. - e2) )**0.5

B1 = a * B * k * (1. - e2)**0.5 / ( 1. - e2 * (SIN(lat0))**2. )

t0 = TAN(pi / 4. - lat0 / 2.) / ( (1. - e * SIN(lat0)) / (1. + e * SIN(lat0) ) ) ** (e/2.)
D = B * (1. - e2)**0.5 / ( COS(lat0) * ( 1. - e2 * (SIN(lat0))**2.)**0.5 )
D2 = D*D
IF (D < 1.) THEN
  !D = 1.
  D2 = 1.
END IF
F = D + (D2 - 1.)**0.5 * SIGN (1.,lat0)
H = F * t0**B
G = (F - 1. / F) / 2.
g0 = ASIN(SIN (azimuth) / D)
l0 = lon0 - (ASIN(G * TAN(g0))) / B
!uc = (B1 / B) * ATAN( (D2 - 1.)**0.5 / COS(azimuth) ) * SIGN (1.,lat0)

!gc = SIN(azimuth) / D
gc = azimuth

t = TAN(pi /4. - lat / 2.) / ( (1. - e * SIN(lat)) / (1. + e * SIN(lat)) )**(e/2.)
Q = H / t**B
S = (Q - 1. / Q) / 2.
TT = (Q + 1. / Q) / 2.
VV = SIN(B * (lon - l0))
UU = (-VV * COS(g0) + S * SIN(g0)) / TT
v = B1 * LOG( (1. - UU) / (1. + UU) ) / (2. * B)
u = B1 * ATAN2( (S * COS(g0) + VV * SIN(g0)) , COS( B * (lon - l0)) ) / B

x = v * COS(gc) + u * SIN(gc) + falseE
y = u * COS(gc) - v * SIN(gc) + falseN

END SUBROUTINE ConvertGeodeticToHotine